Nuclear quantum effects in solids using a colored-noise thermostat 
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We present a method, based on a non-Markovian Langevin equation, to include quantum correc- 
tions to the classical dynamics of ions in a quasi-harmonic system. By properly fitting the correlation 
function of the noise, one can vary the fluctuations in positions and momenta as a function of the 
vibrational frequency, and fit them so as to reproduce the quantum-mechanical behavior, with min- 
imal a priori knowledge of the details of the system. We discuss the application of the thermostat 
to diamond and to ice Ih. We find that results in agreement with path-integral molecular dynamics 
can be obtained using only a fraction of the computational effort. 
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Nuclear quantum effects arc extremely important in 
many condensed-phase systems. For instance, zero-point 
fluctuations affect static correlations, and energy quan- 
tization causes deviations from the classical value of the 
specific heat at low temperatures. A quantum treatment 
of the ionic degrees of freedom is mandatory to capture 
such effects, which are particularly important when light 
atoms or stiff vibrational modes are present. However, in- 
cluding quantum effects is computationally demanding, 
even if one is interested only in equilibrium properties, as 
is the case here. For this reason, the nuclei in molecular 
simulations are often treated classically, even when the 
electronic degrees of freedom are accounted for quantum- 
mechanically [![. 

When exchange-symmetry effects are not relevant, the 
method of choice for studying equilibrium expectation 
values is path-integrals molecular dynamics (PIMD) 0, 
0]. However, this comes at a high computational cost, 
since many replicas of the system need to be simulated 
in parallel. Approximate but less expensive methods such 
as Feynman-Hibbs effective potentials Q have also been 
used, as well as semiclassical approaches to treat zero- 
point energy (ZPE) Their range of validity is limited 
to weak quantum behavior, and to cases where the Hes- 
sian of the potential is available and cheap to compute. 
The interest in methods to introduce quantum effects 
in classical trajectories is thus still very high, see e.g. 
Ref. for recent applications. 

In a recent work 0] we have demonstrated that a gener- 
alized, linear Langevin equation with colored noise can be 
used to obtain a highly tunable thermostat for constant- 
temperature molecular dynamics. In this Letter we will 
show that, by suitably extending this idea, it is possible 
to modify Hamilton's equations so as to introduce the 
quantum-mechanical effects. Furthermore, this comes 
with only a negligible increase in computational effort 
with respect to traditional classical simulations, and re- 
quires a minimal prior knowledge of the properties of the 
system. The idea can be traced back to the semiclassi- 



cal approximation to the quantum Langevin equation [8|] , 
which has been used to model quantum systems in con- 
tact with a quantum harmonic bath. A similar idea has 
been recently used by Buyukdagli et. al [§], in order 
to compute quantum specific heat in harmonic systems. 
However, Buyukdagli's scheme is only qualitatively cor- 
rect even for a harmonic oscillator, and neglects ZPE 
completely. In this Letter we propose a more general 
approach, which allows one to obtain high accuracy in 
reproducing quantum specific heats, and also tackles the 
more challenging task of introducing zero-point motion 
effects. This is achieved by effectively and automatically 
enforcing the quantum position and momentum distri- 
butions. Our method gives excellent results in systems 
with limited anharmonic coupling, and is ideal for treat- 
ing quantum effects in solids. 

Let us first consider a harmonic oscillator which is 
evolved in time according to a generalized Langevin equa- 
tion. This equation can be written in a Markovian form 
by suitably extending the state vector 0, [13, Hi- For 
each degree of freedom, a set of n additional momenta 
Si are introduced, which complement the position q and 
the physical momentum p. The value of n depends on 
the structure of the memory kernel for the generalized 
Langevin equation and, in our experience, a choice be- 
tween 4 and 12 allows sufficient flexibility. For simplic- 
ity, we assume mass-scaled coordinates, q ^ q^/rii and 
p <— p/ y/m. We introduce a compact convention to rep- 
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Thus, a matrix without subscript acts on the subspace of 
additional momenta only. The p and qp subscripts denote 
matrices which also act on the p and on the {q,p) respec- 
tively. The Markovian form for the generalized Langevin 
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equation reads 




where ^ is a vector of n + 1 uncorrelated gaussian random 
numbers. Equation ([2]) has been chosen as the most gen- 
eral Unear stochastic equation where the position q is not 
coupled with the noise nor with the si . This form allows 
for an easier generalization to the anharmonic case. The 
static covariance matrix C^p can be obtained by solving 
the matrix equation [ll|, fl3 | 



AqpCqp + CqpAqp — BqpBgp 



(3) 



In Ref . 



we have chosen Bp so as to obtain Cp = ksT 
{cqq = kBT/uj^, in our units), which corresponds to en- 
forcing detailed balance. In a quantum oscillator at fi- 
nite temperature the distribution of position and momen- 
tum is still Gaussian, but its variance has a nontrivial 
dependence on w, i.e. (p^) = (g^) = ^ coth ^j^j, . 
We can therefore perform a fitting procedure, tuning Bp 
and Ap so that the w-dependence of Cqq — (q^) juP' and 
Cpp = (p^) reproduces closely the exact quantum fluc- 
tuations. The fitting procedure is not trivial, and will 
be discussed elsewhere. Once a set of parameters has 
been found which guarantees a good fit over a certain 
frequency range (wmin, Wmax), the thermostat automati- 
cally enforces the correct quantum fluctuations for any 
system whose typical vibrations fall within this range. 
Tipically, an error below a few percents can be obtained 
over a frequency range spanning several orders of mag- 
nitude. In the harmonic limit, the momentum and posi- 
tion distributions computed using our thermostat sample 
the quantum distributions within the accuracy of the fit. 
Furthermore, it is easy to show that the averages of any 
local operator which depends separately on positions or 
momenta are also correctly computed. 

In the general, anharmonic case, the coupling between 
position and momentum is nonlinear. Therefore, one 
needs to use a finite time step; in particular, one pro- 
ceeds by first integrating the (p, s) variables, which can 
be evolved using the exact finite-time propagator built 
from Ap and Bp, and then updating q and p by integrat- 
ing Hamilton's equations [l3|. 

As a first example, we apply this method to a one- 
dimensional, anharmonic potential of the form 
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For a fixed value of k and temperature T, this poten- 
tial allows one to explore a range of behaviors that goes 
from the highly anharmonic, classical limit at ^ to 
a quantum, harmonic regime at high cu. In Figure [T] we 
compare the exact, quantum solution with the averages 
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FIG. 1: (color online) Mean total ((£)), potential {{V)) and 
kinetic {{K)) energy for a proton in the external potential of 
Eq. © as a function of iv, for A; = 1 and T = 100 K. In 
the inset, the Fourier transform of the momentum distribution 
is reported for the fully quantum, classical and colored-noise 
thermostat simulations at tJ = 200 cm~^. Dashed black lines 
correspond to the exact quantum solution, the dark (red) se- 
ries are results from present work, using a set of parameters 
fitted over a frequency range from 2 to 2000 cm~^. Light 
(blue) lines correspond to the classical expectation values. 



obtained using our colored- noise thermostat. There is 
a remarkably good quantitative agreement not only in 
the asymptotic w — > and a; — > oo limits, but also in 
the intermediate region, where quantum-mechanical and 
anharmonic effects are significant, suggesting that both 
can be captured, albeit not fully. Also the momentum 
distribution (shown in the inset) is in good agreement 
with its quantum-mechanical counterpart. This is par- 
ticularly appealing, since conventional PIMD can only 
sample positions, and one must introduce special proce- 



dures to sample momenta as well |14l . 115 



We now move on to more complex and realistic appli- 
cations. As we have already discussed 0, if one appHes a 
separate thermostat to each degree of freedom, the multi- 
dimensional harmonic problem can still be solved analyti- 
cally, by projecting the dynamics onto the normal modes. 
One would then expect independent phonons to thermal- 
izc at different effective classical temperatures accord- 
ing to the target Cqq {lj) and Cpp (uj) relations, provided 
that the response to the thermostat is faster than the 
anharmonic coupling between different phonons. Heuris- 
tically, one would expect that the error on the energy of 
a phonon of frequency Ui would grow with the coupling 
time T^f(wi), as defined in Ref. [7|, and with the largest 



difference in target energy, AE = Cpp{uj, 



Any such error should instead decrease with the classical 
lifetime of the phonon T^iuji), which gauges the internal 
energy transfer to other vibrational modes. This effect 
can be reduced to a great extent by modifying the fit 
such that not only are the quantum distributions repro- 
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FIG. 2: (color online) Lattice parameter of diamond as a func- 
tion of temperature, computed by isothermal-isobaric runs 
at atmospheric pressure [13]. Black dashed lines correspond 
to PIMD results from Ref. [3], and red circles to averages 
computed from runs using our properly fitted colored-noise 
thermostat. In inset (a) the average kinetic temperature T*, 
projected on a few selected normal modes, is plotted for three 
target temperatures (from top to bottom, 1000 K, 500 K and 
200 K) as a function of the mode frequency. The continuous 
lines are the expected quantum T* (lS) curves. In the inset (b) 
we compare the total internal energy (kinetic plus potential) 
as computed with PIMD (black dashed line, Ref. [3|) and 
with colored-noise thermostat fitted to the quantum E{uj), 
with (red circles) and without (gray lozenges) ZPE. The lat- 
ter (gray) points have been aligned to the PIMD results at 
T ^0. 



duced, but also the correlation time th is made as small 
as possible. 

In order to assess the accuracy of this approach, wc 
first study diamond, an archetypical quasi-harmonic sys- 
tem. We used the Tersoff classical potential [l^, [2^, for 
which accurate PIMD results have been reported [18[ . In 
Figure [2] wc compare observables computed with PIMD 
and with our colored- noise thermostat. The correlated 
Langevin dynamics is able to reproduce quantitatively 
the quantum effects on the thermal expansion and on 
specific heat down to T « 0.16_d. Only at lower temper- 
atures do we start observing significant discrepancies. 

To understand the reason for this breakdown, it is in- 
structive to look at the kinetic temperature T* of the dif- 
ferent phonons [inset (a)]. This is a powerful probe of the 
accuracy of the method, which can be performed when- 
ever a normal-modes analysis is meaningful. In this case, 
the normal-modes analysis shows that, in the extremely 
quantum regime for T < O.lQu, the thermostat fails to 
counterbalance the phonon-phonon coupling due to an- 
harmonicities. Because of this internal coupling, energy 
flows from the stiff modes (which thermalize at a lower 
temperature than expected) to the slow ones (which turn 
out to be hotter than desired). A large energy difference 
AE Ri h (wmax — '^min) /2 must bc maintained between 
fast and slow modes. Moreover, because of ZPE, the 
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FIG. 3: (color online) Radial distribution functions for ice at 
T = 220 K. Black, dashed line corresponds to the PIMD re- 
sults of Hernandez et al. [2l| . dark (red) and light (blue) series 
to colored-noise and purely classical runs, respectively. Runs 
have been performed in isothermal ensemble for a box of 360 
water molecules at experimental density [S^l , and are 125 ps- 
long. In inset (a), we show the experimental proton momen- 
tum distribution for ice at 269 K (black dashed line) [l^ , and 
corresponding distribution as computed from a classical (blue 
line) and colored-noise (red line) simulation using a flexible 
TIP4P-like potential 0| at the rescaled temperature (247 K, 
i.e. 4K below the quantum melting point of the potential). 



motion of the ions is not limited to the harmonic re- 
gion of the potential energy surface. This leads to higher 
anharmonic coupling and eventually to shorter phonon 
lifetimes. This explanation is supported by the results 
shown in inset (b). Here we performed the fit by con- 
sidering an energy-versus-frequcncy relation which ex- 
cludes the zero-point contribution, as done in Ref. 0]. 
The agreement with PIMD energies is now virtually per- 
fect, since high-frequency modes are basically frozen and 
normal modes are almost perfectly decoupled, so that 
internal energy transfer becomes negligible. However, 
if one is interested in quasi-harmonic effects, or simply 
in observables which arc functions of the atomic coor- 
dinates or velocities, neglecting ZPE would mean sam- 
pling an unphysically cold system, where stiff modes are 
completely frozen. Some examples of such observables 
arc the Debye- Waller factor, radial distribution functions 
and momentum distributions. 

To demonstrate the ability of our method to compute 
structural properties, we have performed simulations on 
a TIP4P [ill model of ice Ih at 220 K, which poses the ad- 
ditional challenge of showing larger anharmonicity than 
diamond. In Figure [3] we show the radial distribution 
function goo obtained using a classical simulation and 
our scheme, taking as a reference literature PIMD re- 
sults [21}. The agreement is very good, and the peak 
broadening is correctly and quantitatively predicted. Fit- 
ting the thermostat removing ZPE as in [9|] leads to an 
unphysical reduction of the width of the peaks, which are 
much narrower than for the classical simulation. 
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We also performed similar calculations using a related 
flexible water forcefield [2^, in order to asses the accu- 



racy in a system with an extreme spread of vibrational 
frequencies. Results are in good qualitative agreement 
with rigid- water ones, but we cannot perform a quanti- 
tative comparison in the lack of PIMD results for flexible 
water. This agreement demonstrates that a large AE is 
not a problem, provided that phonon-phonon coupling 
is weak. For this model we also computed the proton- 
momentum distribution (see inset of Figure [3|) . A quan- 
titative comparison with experiment is difficult because 
of the limitations of the model potential, which makes it 
difficult to discuss the origins of the discrepancies. For 
example, the simulation has to be performed 20 K below 
the experimental temperature, because of the dicrepancy 
in the melting temperature. That said, the improvement 
over the purely classical results is impressive, especially 
considering that it has been obtained at negligible com- 
putational cost. 

In conclusion, we have presented a thermostatting 
strategy which can be applied to a vast class of solid- 
state structural problems, including disordered systems 
and glasses. The approach is appealing for several rea- 
sons: no detailed information on the system is required 
(for instance, the same parameters could be used for dif- 
ferent polymorphs of ice), and the correct position and 
momentum distributions are automatically enforced. In 
a semiclassical sense, it can be seen as a method to auto- 
matically equilibrate different normal modes at the ap- 
propriate, frequency-dependent temperature. The im- 
plementation is straightforward, as one only needs to act 
on the velocities, just as with a traditional stochastic 
thermostat. Most importantly, the additional computa- 
tional cost is only noticeable for simulations employing 
simple, short-range, two-body potentials. We arc consid- 
ering strategies to extend the range of applicability to ex- 
tremely anharmonic systems such as liquids. Among the 
possible approaches, the connections between the quan- 
tum Langevin equation and path-integrals formalism sug- 
gest the possibility to combine our method with PIMD. 
Finally, we also notice that virtually any energy-versus- 
frequency curve can be reproduced, so the method can 
be used in other applications beside the simulation of 
quantum effects. 
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